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Abstract 

The multiconfigurational time-dependent Hartree method (MCTDH) [H.-D. Meyer, U. Manthe, 
and L. S. Cederbaum, Chem. Phys. Lett. 165, 73 (1990); U. Manthe, H.-D. Meyer, and L. S. 
Cederbaum, J. Chcm. Phys. 97, 3199 (1992)] is celebrating nowadays entering its third decade of 
tackling numerically-exactly a broad range of correlated multi-dimensional non-equilibrium quan- 
tum dynamical systems. Taking in recent years particles' statistics explicitly into account, within 
the MCTDH for fermions (MCTDHF) and for bosons (MCTDHB), has opened up further op- 
portunities to treat larger systems of interacting identical particles, primarily in laser-atom and 
cold-atom physics. With the increase of experimental capabilities to simultaneously trap mixtures 
of two, three, and possibly even multiple kinds of interacting composite identical particles together, 
we set up the stage in the present work and specify the MCTDH method for such cases. Explicitly, 
the MCTDH method for systems with three kinds of identical particles interacting via all combina- 
tions of two- and three-body forces is presented, and the resulting equations-of-motion are briefly 
discussed. All four possible mixtures (Fcrmi-Fermi-Fermi, Bose-Fermi-Fermi, Bose-Bose-Fermi and 
Bose-Bose-Bose) are presented in a unified manner. Particular attention is paid to represent the 
coefficients' part of the equations-of-motion in a compact recursive form in terms of one-body den- 
sity operators only. The recursion utilizes the recently proposed Combinadic-based mapping for 
fermionic and bosonic operators in Fock space [A. I. Strcltsov, O. E. Alon, and L. S. Cederbaum, 
Phys. Rev. A 81, 022124 (2010)] and successfully applied and implemented within MCTDHB. Our 
work sheds new light on the representation of the coefficients' part in MCTDHF and MCTDHB 
without resorting to the matrix elements of the many-body Hamiltonian with respect to the time- 
dependent configurations. It suggests a recipe for efficient implementation of the schemes derived 
here for mixtures which is suitable for parallelization. 

PACS numbers: 31.15.xv, 67.60.-g, 05.30.Fk, 05.30.Jp, 03.65.-w 
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I. INTRODUCTION 



Quantum non-equilibrium dynamics is important to many branches of physics and 
chemistry [1-6] and often requires the solution of the time-dependent many-particle 
Schrodinger equation. A particular efficient method to solve the time-dependent many- 
particle Schrodinger equation is the multiconfigurational time-dependent Hartree (MCTDH) 
algorithm and approach [7-10]. MCTDH, which is considered at present the most efficient 
wave-packet propagation tool, has amply been employed for multi-dimensional dynami- 
cal systems of distinguishable degrees-of-freedom, typically molecular vibrations, see, e.g., 
Refs. [11-18]. We mention that recent developments on multi- layer formulation of MCTDH 
have opened up further possibihties to treat larger systems of distinguishable degrees-of- 
freedom [19-21]. MCTDH has recently been apphed with much success to various systems 
with a few identical particles in the field of cold-atom physics, see, e.g., Refs. [22-27]. 

In recent years, taking the quantum statistics between identical particles a priori into 
account, the MCTDH method has been specified for systems of identical particles, which 
opened up interesting possibilities to treat larger systems. First MCTDHF - the fermionic 
version of MCTDH - was developed by three independent groups [28-30]. Shortly after, 
MCTDHB - the bosonic version of MCTDH - was developed in [31, 32]. For applications of 
MCTDHF to laser-matter interaction and other few-fermion problems see, e.g., Refs. [33- 
42], where the last work combines optimal control theory with MCTDHF. For applications 
of MCTDHB to Bose- Einstein condensates see, e.g., Refs. [43-47], where the last two works 
combine optimal control theory with MCTDHB. 

Since the seminal paper of Lowdin [48], reduced density matrices and particularly re- 
duced two-body density matrices have been a lively field of research, sec, e.g., Refs. [49-55]. 
Reduced one-body density matrices are an inherent part of the MCTDH [7-10]. In the 
present context, reduced one- and two-body density matrices were first used to derive the 
static self-consistent theory for bosons, the multiconfigurational Hartree for bosons (MCHB) 
in [56]. Thereafter, MCTDHB and MCTDHF were formulated in a unified manner by em- 
ploying reduced one-, two- [57] and three-body [10] density matrices. Further specification 
of MCTDH to mixtures of two kinds of identical particles (MCTDH-FF for Fermi-Fermi 
mixtures; MCTDH-BF for Bose-Fermi mixtures; and MCTDH-BB for Bose-Bose mixtures) 
was put forward in [58]. All the above developments made use of the fact that the mean- 



3 



field operators in the traditional MCTDH can be factorized to products of reduced density 
matrices times one-body operators. Finally, we mention that MCTDH has been extended to 
systems with particle conversion (termed MCTDH- conversion), where particles of one kind 
can convert to another kind [59]. 

A breakthrough in the formulation [60, 61] and implementation [62] of MCTDHB has 
stemmed from a general Combinadic-based mapping of bosonic (and fermionic) operators 
in Fock space. In this formulation, the direct calculation of the matrix representation of 
the Hamiltonian in the (huge) multiconfigurational space is abandoned, and is replaced by 
the action of one-body and two-body density operators on the multiconfigurational wave- 
function. The operation of the various density operators can be performed in parallel [62], 
which further accelerates the performance of the algorithm. This brings us closer to the 
topic and contents of the present work. 

Two-body interaction is the most basic interaction in an interacting (quantum) system. 
When the particles comprising the quantum system have internal structure, higher-order 
interactions (forces) may come into play. For instance, in nuclear physics it has long been 
accepted that three-body interactions are necessary to fully understand the structure of 
nuclei, see, e.g. [63, 64]. Much more recently, and in the context of another field, the 
proposition to utilize cold polar molecules to engineer (condensed-matter) systems with 
three-body interactions has been made [65]. So, the motivation to study the non-equilibrium 
dynamics of systems with up to three-body forces is clear. 

But why study the quantum dynamics of a mixture of three kinds of identical particles? 
Are such systems present in nature? In the cold- atom world, the plurality of atoms is one of 
the most important ingredients experimentalists (and theorists) have at their disposal. For 
instance, the element Yb has seven stable isotopes (5 bosonic and 2 fermionic isotopes). Yb 
has been envisaged to play an instrumental role in realizing various interesting ultra-cold 
mixtures (see Ref. [66] for a realization of a Bose-Einstein condensate with ^™Yb atoms 
and the discussion therein). More recently, a quantum degenerate Fermi- Fermi mixture of 
^Li-^°K atoms coexisting with a Bose-Einstein Condensate of ^^Rb atoms were reahzed [67], 
as well as a triply quantum-degenerate mixture of bosonic "'^K atoms and two fermionic '^''K 
and ^Li atoms [68]. Hence, mixtures of three kinds of identical particles have been created 
in the lab. 

All the above dictate the purposes and contents of the present work. The MCTDH 
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method for systems with three kinds of identical particles interacting via all combinations 
of two- and three-body forces is derived, and the resulting equations-of-motion are briefly 
discussed. All four possible mixtures (Fermi- Fermi- Fermi, Bose-Fermi-Fermi, Bose-Bose- 
Fermi and Bose-Bose-Bose) are presented in a unified manner. Particular attention is paid 
to representing the coefficients' part of the equations-of-motion in a compact recursive form 
in terms of one-body density operators only. The recursion utilizes the recently proposed 
Combinadic-based mapping [60] which has already been successfully applied and imple- 
mented within MCTDHB [62]. Our work sheds new light on the representation of the 
coefficients' part in MCTDHF and MCTDHB without resorting to the matrix elements of 
the many-body Hamiltonian with respect to the timc-dcpcndcnt configurations, and sug- 
gests a recipe for efficient implementation of the theory derived here for mixtures which is 
suitable for parallelization. 

The structure of the paper is as follows. In Sec. II we present the building bricks of 
the theory by reconstructing MCTDHF and MCTDHB. In Sec. HI we assemble from these 
ingredients the multiconfigurational time-dependent Hartree method for mixtures of three 
kinds of identical particles interacting via up to three-body forces. A brief summary and 
outlook are given in Sec. IV. Finally, we collect in Appendixes A-C for completeness and ease 
of presentation of the main text various quantities appearing and needed in the derivation. 
The paper and the Appendixes are detailed and intended also to serve as a guide for the 
implementation of the equations-of-motion. The reconstruction of MCTDHF and MCTDHB 
is given in sufficient detail. This allows us to defer to the Appendixes much of the lengthly 
formulas used later on for the mixtures. 
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II. BUILDING BRICKS: RECONSTRUCTING MCTDHF AND MCTDHB 



A. From basic ingredients to mapping 

Our starting point is the many-body Hamiltonian of Na interacting identical particles of 
type A: 

+ IJd^' *^(x)*^(x')H^(^)(x,x')*^(x)*a(x) + (1) 

where h^^^ is the one-body part, W^(^) the two-body part and U^^^ the three-body part. The 
operators h^'^\ W^(^) and U^"^^ can generally be time-dependent. 

We use the time-independent field operator expanded by time-dependent orbitals: 

*^(x) = ^afc(t)0fe(x,t), (2) 

k 

where the annihilation and creation operators obey the usual fermionic/bosonic 
anti/commutation relations, aq(t)a|,(t)±a|,(t)aq(t) = 6kq. Correspondingly, the field operator 
obeys the anti/commutation relations, ^a(x) |^a(x')| ± |^^(x')| ^a(x) = 5(x — x'). 
Here and hereafter the upper sign refers to fermions and the lower to bosons. The coor- 
dinate X = {r,(j} stands for spatial degrees of freedom and spin, if present. Thus, the 
shorthand notations 5(x — x') = 5{r — r')5cr,cr' and f dx = J drY^^ are implied throughout 
this work. Furthermore, we do not denote explicitly the dependence of quantities on time 
when unambiguous. 

Plugging the expansion (2) into the many-body Hamiltonian (1) one gets: 

^{A) ^ ^ + - ^ ^LqlPkslg + g ^ ^kslqlrPkslrlq^ (3) 

k,q k,s,q,l k,s,p,r,l,q 

where the matrix elements with respect to the orbitals {0jk(x, t)} are given by: 
C = / </'I.(x,t)/i(^)(x)0,(x,t)dx, 

= j j </>^(x,t)0:(x',t)W^(^)(x,x')0,(x,i)c/.,(x',t)rfxdx', 

<U - /// 0:(x,t)^:(x',t)0;(x",i)t/('^)(x,x',x")x 

X (j)g{x,t)(f)i{x',t)(j)r{:x.",t)dxdx'dx". (4) 



In (3), we introduce the one-body density operators 



Pkq = 4«9' (5) 

as well as the two- and three-body density operators 

Pkslq - O'kO^s'^l'^q - ^Pkg ^sl =F Pkl Psq 

Pksprlq - dk^^s^^p^rdiaq - ^Pkslq^^pr " Pksrq^l + PksrlPpq ' l^j 

The reason for this choice of notation with density operators in (3) will become clear below. 

4) ~ 
slq 



Wc sec that the two-body density operators can be written as products of the one- 

body density operators, and that the three-body density operators |/0^'^prig| ^® written 
as products of the two- and one-body density operators, and so on, recursively. Hence, the 
one-body density operators in (5) are our basic building bricks. 

The many-body wave-function is expanded by time-dependent configurations (determi- 
nants \i\t) for fermions, permanents |n;t) for bosons) assembled by distributing the Na 
particles over the Ma time-dependent orbitals introduced in the expansion (2). For fermions 
we write [60] : 

con J 

{i} ^A=l 

where the address J a is defined as follows: 

whereas for bosons we write [60]: 

conj 

|*(^)(t)> = J2C4t) \n;t) ^ ^^-W l-^^;^) ' (9) 

{n} Ja=1 



where the address Ja is defined as follows 

Ma-1 



/../.(n) = l+ Y. (10) 

The notation used in (7-10) follows the Combinadic-based addressing scheme of config- 
urations introduced in [60]. For fermions we enumerate configurations by holes, i = 
{ii,i2,---,ij = q, . . . ,iMA-NA) and i''* = {ii,i2, . . . ,ii = k, . . . ,iMA-NA), whereas for 



bosons we enumerate configurations by particles, n = {rii, . . . ,nk, ■ ■ ■ ,nq, . . . , uma) 
n'^i — (m, . . . , rife — 1,. . .,nq + 1, . . . , hma)- The index J a is termed "address" because 
it is an integer uniquely identifying a configuration which is described by the positions of 
the holes i (for fermions) or the occupation numbers n (for bosons). For more details of the 
Combinadic-based mapping and particularly the connection between the bosonic occupation 
numbers and the positions of the fermionic holes see [60]. 

For our requirements, we will need the result of the operation of the basic building bricks 
onto the state vector, namely, the operation of the one-body density operators onto 
|*(^)(t)>. Thus we have: 



conj conj / a \ 



Ja=1 Ja=1 



For fermions we have the following relations 



Cj.^ijit); k^q,k^i , (12) 

0; otherwise 

where the distance between the ij-th hole of i at orbital q and the irth hole of i'^^ at orbital 
k is given by d{i'"^) — \k — q\ — \j — l\ — l [equivalently, d{i'''^) — J2pe{k q) '^v simply enumerates 
how many fermions are there between the A;-th and g-th orbitals]. For bosons we have the 
following relations [60]: 



.(A) ,(A) 

c7: (t) - c'jX)it) = \ ' """'"^"^ ■ v-^ V ■ - — , (13) 

CjA{n){t) xuk; k^q 



CjAinki){t) X y/nj^y^Uq + 1; k^q 



which concludes our exposition of the Combinadic-based mapping and assembly of the op- 
erations of the basic building bricks |pfc^^| on the many-body wave-function |^'*^"^)(t)). 
From Eqs. (5,6) we see how to use the one-body (basic) building bricks to assemble 

higher-body operators. In particular we find: 

c;r"'(o = ^^vrC'jTit) - Spic'/rii) + (0- (14) 

The meaning of the two levels of density operators in the superscripts of the coefficients 

.(A) .(A) 

Cj*'''(t) and Cj*''"'''(i), resulting from higher-body operators in (14), is that the lower-level 



density operator is multiplied on the many-body wave-function first, and the upper-level 
density operator is multiplied thereafter on the result. 

The key ingredient in the utilization of the Lagrangian formulation [32, 69, 70] of the 
(Dirac-Frenkel [71, 72]) time-dependent variational principle to derive the equations-of- 
motion is the evaluation of matrix elements with respect to the multiconfigurational wave- 
function This will be utilized in the next subsection JIB. For the moment, 
we would like to prescribe how such matrix elements with respect to are to be 
evaluated. 

Consider the operator 0*^"^^ which can be a one-body operator, two-body operator, three- 
body operator, etc. Then, we express and compute the expectation value of O^^^ with respect 
to \^^^\t)) as follows [60]: 



conj 



Ja=1 



where 



TV 



(A) 

conf 



N 



(A) 
conf 



Ja=1 



Ja=1 
(A) 4 A) 



In particular, for a one-body operator, O^^^ = Ylik q^kq Pkq ■> §^^- 



(16) 



Ma 



k,q 

for a two-body operator, 6(^) = \ Y.k,s,q,i Oi%pi%: we get from (14): 



(17) 



1 (A) 



k,s,q,l 
T Ma 

k,s,q,l 



(A) 



'(A) 



^SsiC;: it) T c'j7 it) 



(18) 



and for a three-body operator, O(^) = \ T.k,s,p,rM^kspqiAspriq^ we get from (14): 

Ma 

6 

Ma 



1 (A) 

^Ja i^) - a '-'kspqlr^JA 



it) 



k,s,p,r,l,q 



g / J kspqlr 
k,s,p,r,l,q 



(19) 



Finally and generally, the result of a sum of (operations of) operators, e.g., O^^^ + 02^\ on 
|^'(^)(t)) translates to the sum of the respective coefficients [60]: 

A(-4) , A(A) A(A) A(A) 

CjI ^""^ {t)^C^l {t) + C^^ (t). (20) 

These compact relations resting on one-body density operators only [the two-body density 
operators in (19) are assembled from one-body density operators according to (5,6)] will 
be used to reformulate MCTDHF and MCTDHB in a recursive manner in the following 
subsection II B. 



B. Equations-of-motion utilizing one-body density operators and Combinadic- 
based mapping 



We can derive (reconstruct) the MCTDHF and MCTDHB equations-of-motion, taking 
into account a-priori that matrix elements of the form of (15) enter the variational formu- 
lation. Within the Lagrangian formulation [32, 69, 70] of the (Dirac-Frenkel [71, 72]) time- 
dependent variational principle, the action functional of the time-dependent many-particle 
Schrodinger equation takes on the following form: 

d 



s[{Cj,{t)},{M^,t)}] 



dt\ (^^^^\t) 



dt 



Ma 



k,j 



N 



(A) 

conf 



Ja=i 



(21) 



where the timc-dcpcndcnt Lagrange multipliers ^'^'c introduced to guarantee the 

orthonormality of the orbitals at all times. Furthermore, they enable one to first evaluate the 
expectation value of H^^^ — with respect to \^^^\t)) and then to perform the variation, 
which is precisely what is needed in order to exploit the Combinadic-based mapping [60] 
a-priori in the derivation of the equations-of-motion. The (here redundant) time-dependent 
Lagrange multiplier e^^^ (t) ensures normalization of the expansion coefficients at all times, 
and would resurface in the static theory in the case of the imaginary-time-propagation 
formulation. 

To perform the variation of the action functional with respect to the coefficients, we 



express the expectation value ( ^^^^ (t) 



'dt 



^^^\t) ) following the Combinadic-based 
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mapping [60] and the compact expression in Eq. (15): 



AT' 



(-4) 



Ja=1 



C 



Ja 



: a (A) 



(22) 



Representation (22) makes it clear what the variation with respect to the coefficients 
{C}^{t)] would lead to. When this variation is performed exphcitly, one immediately finds: 



C 



Ja 



a_(A) 
at 



(23) 



The meaning of i-§i^^^ is that the time-derivative is a one-body operator in the A-species Fock 
(and orbital) space. According to the rules of the previous subsection II A, the left-hand-side 
of Eq. (23) is given by the sum of its one-, two- and three-body constituents: 



cl it) 



c 



(24) 



The invariance of |\l/*^^^(t)) to unitary transformations of the orbitals, com- 
pensated by the 'reverse' transformations of the orbitals is well-known [7, 8, 
32] and can be represented as follows: 1*^^^^)) = ^jlZiCj^it) \JA;t) = 



Cj^(t)| Ja; t), with obvious notation. This invariance can be utihzed to bring the 
equations-of-motion into a simpler form (see, in particular, the discussion below on the or- 
bitals' part). Primarily, the differential conditions first introduced by the MCTDH founders 

[7, 8]: 



(f>k 



(25) 



kq 



come out explicitly from such a unitary transformation [32, 59] and straightforwardly lead 
in the case of the equations-of-motion for the coefficients to: 



cff{t) = iCjAt). 



'Ja 
'Ja 



cr{t) + cy-'{t) + cz'{t). 



(26) 



For the general form of the differential conditions, Eq. (25), see the literature [9, 10]. We 
remark that a particular interesting representation (put forward and utilized so far for dis- 
tinguishable degrees-of-freedom only) of the differential conditions can be made in order to 
propagate the systems' natural orbitals [73, 74]. 

In MCTDHF and MCTDHB the integration of the coefficients' part in time is performed 
(for unitary time-evolution) by the short iterative Lanczos (SIL) algorithm [75]. We remark 
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on the numerical implementation of Eq. (26) within SIL propagation [62] . For the SIL one 
needs to operate with powers of H onto the many-particle wave-function and construct 
the X-dimensional Krylov subspace: | I^^^H*)) , ^^^^ |*^^^(^)> H^^^^ ^ |^(^)(f)> |. 
In the language of the Combinadic-based mapping of coefficients and utilizing the recipe 
of how to operate with operators on the many-particle wave-function discussed above [60], 



this construction translates to: •|Cj^(t), C_^'^' (t), C^^^ (t), ■ ■ ■ 

Let us now move to the equations-of-motion for the orbitals {(j)k{:x.,t)}. For this, the 
expectation value of the many-body Hamiltonian If^^^ with respect to has to be 

expressed in a form which allows for variation with respect to the orbitals, namely as an 
explicit function of the quantities (integrals) h^j^\ W^f^i and U^f^^j^ in (4). The result reads: 



di 



Ma 
k,q=l 



kq 



+ 



(27) 



T Ma 
+ 9 Pk 



(A) T^(A) 
slq ksql 



Ma 



N' 



(A) 

conf 



+ Q ^ Pkslrlq^Llqlr ^ (^) ^-^A (^) • 

k,s,l,q=l k,s,p,r,l,q=l Ja=1 

The expectation values of the density operators p\^^ , p^^j^ and p^^l^ig with respect to | (t) ) 
(resulting from the expectation value of the Hamiltonian with respect to many-particle wave- 
function) are computed following Eq. (15): 



(A) 
Pkq 



conf 

Ja=1 



N 



(.A) 

conf 



N 



(-4) 
conf 



Ja=1 



pit=T.CjAit)Cjrit), 

Ja=1 



,-SA) 



(28) 



where the coefficients Cj^' (t), Cj^'''(t) and Cj^^p-^'H) are given in Eqs. (12,13) and (14), 
respectively. We collect the expectation values of the one-body density operators as the 
matrix p(^)(t) = {pi^)(t)}. 

One should remember that the expectation values of two- and three-body density opera- 
tors can generally not be factorized into products of expectation values of one-body density 
operators. For instance (and in the language of the Combinadic-based mapping of coef- 

(A) 

ficients), c'/;'''{t) = ±SsiC'/; it) T c'/; (t) ^ ±SsiC'/; (t) ^ c'/J m'/; (t). This is 

unlike the operation of the density operators themselves on the many-particle wave-function 
utilized above. 
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We can now perform the variation of S [{C'j^(i)} , {4>k{^, t)}] with respect to the orbitals. 
This variation has been detailed in the hterature, see [32, 57], and we give here the main 
steps in the derivation of the equations-of-motion as far as they are needed for our needs 
later on. Making use of the orthonormality relation between the time- dependent orbitals 
{0fc(x, t)}, we can solve for the Lagrange multipliers, k, j = 1, . . . , Ma- 



(29) 



Ma 



(A) 

q 



d_ 
di 



(A) 



Ma ^ Ma 

Pkslq^sl + 2 ^ Pksprlq^splr 
s,p,r,l=l 



s,l=l 



The Lagrange multipliers |A*i^^(0| can be eliminated from the equations-of-motion which 
is achieved by the introduction of the projection operator: 



Ma 



(30) 



u=l 



When this is done, we find the following equations-of-motion for the orbitals {0j(x, i)}. 



Ma 



(31) 



Ma Ma ( Ma \ 

+ E E ( o'Cl,"^^' + \ E "it-.C ) I*. 

k,q=\ 



s,l=\ 



p,r=\ 



where 



Wf{:^,t) = I 0:(x',t)#(^)(x,x')0Kx',t)<ix', (32) 

u^splM^ 1 1 0:(x',f)0;(x",i)t/(^)(x,x',x")<^Kx',t)0.(x",t)c/x'dx", 

are local (for spin-independent interactions), time-dependent one-body potentials, and 4>j = 



dt 



Utihzing the differential conditions (25) we can ehminate the projection operator P^^^ 
appearing on the left-hand-side of Eq. (31) and arrive at the final result for the equations- 
of-motion of the orbitals in MCTDHF and MCTDHB (see [10, 57]), j = 1, . . . , M^: 



Ma 



(33) 



Ma Ma ( -.Ma \ 

E {p'-mt E f^^^' + 5 E i^lZ:/'^ I* 



k,q=l 



s,l=l 



p,r=l 
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Summarizing, the coupled sets of equations-of-motion (26) for the expansion coefficients and 
(33) for the orbitals constitute the MCTDHF and MCTDHB methods, where the one-body 
density operators (5,6) are employed as the basic building bricks in their construction and 
implementation. 

We can also propagate the MCTDHF and MCTDHB equations-of-motion (26,33) in 
imaginary time and arrive for time-independent Hamiltonians at the corresponding self- 
consistent static theories, MCHF [76, 77] and MCHB [56]. Thus, setting t -it into the 
coupled sets (23,31) or into (26,33), and translating back from the projection operator P^^^ 
to the Lagrange multiphers jA^if^j) ^^e final result reads, k — 1, . . . , Ma' 



Ma 

E 

q=l 



Ma / ^ Ma 

Pkq ^'^^^ + ^ ( Pkslq^sl + 2 ^ Pksprlq^splr 
s,l=l \ p,r=l 



Ma 



Cr = ^^^^Cj,^ ^Ja, (34) 



where, making use of the normalization of the many-particle wave-function, — 
X]j^=i ^Ja^Ja eigen-energy of the system. Making use of the fact that the ma- 

trix of Lagrange multipliers {/i^^'*} is Hermitian (for stationary states) and of the invariance 
property of the multiconfigurational wave-function (to unitary transformations of the or- 
bitals compensated by the 'reverse' transformations of the coefficients), one can transform 
Eq. (34) to a representation where {At^^^} is a diagonal matrix. 

All in all, we have formulated in the present section the MCTDHF and MCTDHB 
equations-of-motion, as well as their static variants MCHF and MCHB, by (i) utilizing 
in a recursive manner one-body density operators only, and by (ii) employing a priori the 
Combinadic-based mapping formulation of Ref. [60] to evaluate matrix elements. This sets 
up the tools to put forward the MCTDH theory for mixtures of three kinds of identical 
particles in the following Sec. HI, and to briefly discuss its structure and properties, and 
how to implement it. 



III. THREE KINDS OF IDENTICAL PARTICLES: MCTDH-FFF, MCTDH-BFF, 
MCTDH-BBF AND MCTDH-BBB 

In the present section we specify the MCTDH theory for mixtures of three kinds of iden- 
tical particles, interacting with up to three-body forces. Before we get into the details of 



14 



derivation and flood of equations, we would like to lay out a general scheme or flowchart 
that one can follow to handle similar or even more complex mixtures. Speciflcally, we need 
to assign a different set of time-dependent orthonormal orbitals to each and every species 
in the mixture. These orbitals are then used to assemble the time- dependent configura- 
tions (with determinants' parts for fermions and permanents' parts for bosons). The many- 
particle wave-function is thereafter assembled as a linear combination of all time-dependent 
configurations with time-dependent expansion coefficients. The many-particle Hamiltonian 
contains different terms: It contains intra-species terms and inter-species terms which con- 
sist of two-body, three-body and so on interactions. The main point in the representation of 
the Hamiltonian is the utilization of one-body density operators. In turn, all intra-species 
and inter-species interactions can be represented utilizing (products of) one-body density 
operators only. 

The key step in the derivation of the equations-of-motion is the utilization of the La- 
grangian formulation [32, 69, 70] of the (Dirac-Prenkel [71, 72]) time-dependent variational 
principle with Lagrange multipliers for each species' orbitals, ensuring thereby the orthonor- 
mality of the orbitals for all times. In such a way, matrix-elements appear within the formu- 
lation explicitly, before the variation with respect to either the expansion coefficients or the 
orbitals is performed. The equations-of-motion for the expansion coefficients of the multi- 
configurational wave-function are obtained by taking the variation of the action functional 
when it is expressed explicitly in terms of the expansion coefficients. The Combinadic-based 
mapping [60] lifts the necessity to work with the huge matrix representation of the Hamil- 
tonian with respect to the configurations, and allows one to efficiently perform operations 
on the vector of expansion coefficients directly. The equations-of-motion for the orbitals 
are obtained by taking the variation of the action functional when it is expressed explicitly 
in terms of the (integrals of the) orbitals. When this is performed, expectation values of 
the various density operators in the Hamiltonian (with respect to the many-particle wave- 
function) emerge which can be efficiently computed utilizing the Combinadic-based mapping 
[60]. 
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A. Additional ingredients for mixtures 



For a mixture of three kinds of identical particles, Na particles of type A, Nb particles 
of type B and Nq particles of type C, we need now two additional field operators expanded 
by different complete sets of time-dependent orbitals: 

*B(y) ^J2^h'{t)iPk'iy,t), *c(z) = (35) 

k' k" 

where the field operator for the A-species particles (x) was first introduced and expanded 
in (2). Note that each species can have a different spin, hence the explicit three distinct 
coordinates x, y and z. Field operators of distinct particles (can be chosen to) commute. 

Our starting point is the many-body Hamiltonian of the most general 3-species mixture 
with up to 3-body interactions: 

^fj{AAB) ^ fj{ABB) ^ fj{AAC) ^ fj{ACC) ^ (j{BBC) ^ (j{BCC) ^ (j{ABC) _ 

Here, H^'^\ H^^^ and H^'-^^ are the single-species Hamiltonians that can be read of (1). The 
inter-species two-body interaction parts are given by: 

W^^^^ = 1 1 dxdy*^(x)*^(y)ll^(^^)(x,y)*B(y)*^(x), 
W^^c) ^ J J rfxdz#^(x)*^(z)iy(^'^)(x,z)*c(z)*A(x), 

^ // ^y^z*L(y)*L(z)^^''''Hy,z)*c(z)*B(y). (37) 

The inter-species three-body interaction parts, resulting from the force between two identical 
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particles and a third distinct particle, are given by: 




fj{ABB) ^ 



(j{AAC) ^ ^ 




(j{ACC) 



(j{BBC) 



(j{BCC) 



1 

2 

x*B(y)*A(x')*A(x), 

1 

2 

x*B(y')*s(y)*A(x), 
1 

2 

x*c(z)*a(x')*a(x), 

1 




dx«!ydy'*^(x)*l,(y)*l,(y')m^^^) (x, y, y') x 



(/x(/x'dz*^(x)*i(x')*5;(z)?7(^^^)(x,x',z) X 



dyidzdz'^\(x.)%(z)%(z')U^^^^\yi, z, z') x 



x*c(z')*c(z)*a(x 

1 




clyc^y'(iz*^(y)*^(y')*J;(z)t/(^^^) (y, y', z) x 



x*c(z)*s(y')*s(y), 

1 

2 

X*c(z')*c(z)*B(y)- 




dydzdz'*Uy)*L(z)*L(z')^^''''^ny,z,z') 



X 



(38) 



Finally, the inter-species three-body interaction part, resulting from the force between three 
different particles is given by: 

^ // ld^dydz^\{^)^l{y)^l{z)m^^''\^,y,z) x 

x*c(z)*B(y)*A(x). (39) 

When all the above are combined, i.e., the field operators ^^(x), ^^(y) and ^b(z) sub- 
stituted into the various interaction terms, we find the following second-quantized expression 
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for the mixture's Hamiltonian: 

~ Z_^"'kq Pkq ^ ^ 2-^ ksqlPkslq ^ q /-^ kspqlryksprlq ^ 

k,q k,s,q,l k,s,p,r,l,q 

^ '^k'q'Pk'q' ^ ^ 2-^ k's'q'l'Pk's'l'q' ^ g ^ k's'p'q'l'r' Pk's'p'r'l'q' ^ 

k',q' k',s',q',l' k' ,s' ,p' ,r' ,1' ,q' 

^ "'k"q"Pk"q" ^ ^ 2-^ k"s"q"l"Pk"s"l"q" ^ 

k",q" k",s",q",l" 



^ k" ,s" ,p" ,r" ,1" ,q" 



^ ^^kk'qq'Pkq Pk'q' ^ ^^kk"qq"Pkq Pk"q" ^ 2-^ k'k"q'q"Pk'q'Pk"q" ^ 

k,k',q,q' k,k",q,q" k',k",q',q" 

+ i V ui^^B) JA) JB) I MABB) JA)JB) 

2-^ '^kk'sqq'lPkslqPk'q' ^ cy 2-^ kk' s'qq'V Pkq Pk' s'l'q' ^ 



2 ^ / KK' sqq' ir KSiqi K' q' ' 2 

k,k' ,s,q,q' ,1 k,k' ,s' ,q,q' ,1' 



^2 '^kk"sqq"lPkslqPk"q" ^ 2 2-^ kk"s"qq"l"Pkq Pk"s"l"q" ^ 

k,k",s,q,q",l k,k" ,s" ,q,q" ,1" 

1 tAbbc) ^b) 4c) 1 tAbcc) 4b) 4c) 

"•"9 ^k'k"s'q'q"l'Pk's'l'q'Pk"q" ' r, 2-^ ^ k'k" s"q'q"l" Pk'q' Pk" s"l"q" 



2 / J K' K" s' q' q" l' r K' s' I' q' r n" q" ' 2 

k' ,k" ,s' ,q' ,q" ,1' k',k",s",q',q",l" 

^ '^kk'k"qq'q"Pkq Pk'q'Pk"q"- \^^) 

k,k',k",q,q',q" 

1{(ABC) governs the non-equihbrium dynamics (and statics) of the mixture, and the most 
efficient way to treat this dynamics is by specifying the MCTDH method for the mixture, 
making use of the building bricks of the previous section II. We see in (40) two kinds of 
ingredients. First, there are matrix elements (integrals) of the various interaction terms 
with respect to the orbitals. For the flow of exposition and for completeness, we list them 
in Appendix C. Second, there are various density operators in H(^^^\ The B and C 
intra-species density operators can be read directly from Eqs. (5,6), when replacing therein 
the A-species quantities. The inter-species density operators in (40) can all be represented 
as appropriate products of the one-body density operators: IpIy} {'^l;"?"}- 

These are the (basic) building bricks of our theory for mixtures. But how to operate with 
them on many-particle wave-functions of mixtures? 

The multiconfigurational ansatz for a mixture of three kinds of identical particles now 
takes on the from: 

jv(^) jv(^) jv(^' 

conj conj conj 

|^(ABC)(^)^ =Y^Y.Y. Cj,,j,,j,{t)\JA,JB.Jc;t), (41) 
Ja=1 ^b=1 ^c=1 
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where we denote hereafter J — {Ja,Jb,Jc) for brevity, such that Cj{t) = Cj^^jg^j^{t), 

J; t) ^ \Ja. Jb, Jc; t) and = Ejl^i ^J^S ^Jc=v 

To prescribe the action of operators on the multiconfigurational wave-function of the 
mixture (41), all we need to know is how the density operators operate on |\l'*^"^^'-^)(t)). The 
operation of the basic, one-body density operators, whether p\^\ p^y or can be read 
of directly from Eqs. (11-19) and we will not repeat them here (one needs just to replace 
therein by J in the overall notation, and Ma by Mb or Mc, when appropriate; also see 
[60]). For the inter-species two-body density operators we have: 

~(A) .(B) .(A) .(C) .(S) .(C) 

C5,'=^''='^'(i), C)''''"'^"{t), cy {t). (42) 

The notation in (42) is to be understood as follows: The two one-body density operators 
(in each case) are written as superscripts on the same level, signifying that they commute 
one with the other; The operation of the two one-body density operators on \'^'^'^^^\t)) is 
to be performed sequentially, i.e., the first operates on |^('^^^)(t)) and the second operates 
on the outcome. Finally, for the inter-species three-body density operators we have: 

.(A) .(B) .(A) .(B) .(A) .(C) 

QpkslqPyq' Q^kq Pk's'l'q'^^^ QpkslqPy'q" 

.(A) .(C) .(B) .(C) .(B) .(C) 

Qpkq Pk"s"l"q" Ij.-^ Qpk's'l'q'Pk"q" ^^-j Qp k' q' P k" s" I" q" ^^-^ ^^g^ 

where the operation of the two-body density operators appearing in the superscripts is 
further decomposed to operations of one-body density operators on \^^'^^^\t)) analogously 
to Eq. (14) [see Appendix A], and 

.(A) .(B) .(C) 

(jPkq Pk'q'Py'q" (44) 

Now we are in the position to write the action of operators on the multiconfigurational 
wave- function of the mixture (41). This is collected for ease of reading and for completeness 
in Appendix A. 

We have gathered most ingredients for the derivation of the cquations-of-motion, which 
is written down in the subsequence section IIIB. There are four possible mixtures (Fermi- 
Fermi- Fermi, Bose- Fermi- Fermi, Bose-Bose- Fermi and Bose-Bose-Bose), and the resulting 
MCTDH-FFF, MCTDH-BFF, MCTDH-BBF and MCTDH-BBB are to be derived and pre- 
sented in a unified manner, in the spirit it has been done in the single-species case [10, 57] 
(and the previous section II) and for mixtures of two kinds of identical particles [10, 58]. 
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B. Equations-of-motion utilizing one-body density operators and Combinadic- 
based mapping for mixtures 



The action functional of the time-dependent many-particle Schrodinger equation takes 
on the form: 



s [{Cjit)} , {0fe(x, t)} , {i^k'iy, t)} , {xk"{z, m = 



(45) 



dt\ 



H(ABC)_-d_ ^iABC)^A_ 



Ma 



dt 

Mb 



k,j k',j' 
Mc 



k",j" 



{f} 



where the time-dependent Lagrange multipliers Ia^I^^I^)!? 

introduced, respectively, to ensure the orthonormality of the A-, B- 
and C-species orbitals at all times. Note that orbitals of distinct particles need not be 
orthogonal to each other. As for the single-species theory, the Lagrange multiplier e^"^^*^) (t) 
is redundant in the time-dependent case and will resurface in the static theory. 

In what follows we present the main steps of the derivation. More details and vari- 
ous quantities needed for the derivation and in particular for the implementation of the 
equations-of-motion are deferred to Appendix B and Appendix C. 

To perform the variation of the action functional (45) with respect to the coefficients, we 
write the expectation value of H^^^^") with respect to \^!^^^^\t)) in a form which is explicit 
with respect to the coefficients: 



at 



C5 (^t)-iCf{t) 



(46) 



The three time-derivative operators i-§i^^\ '^'§1^^^ ^^^^^ make it clear that to each species 
there is associated a different one-body operator representing the derivative of orbitals in 
time. 

Performing the variation of S [{Cj{t)} ,{(f)k{x,t)} ,{i/jk'{y,t)} ,{xk"{'^,t)}] with respect 
to the expansion coefficients |Cj(t)|, we then make use of the differential conditions for 
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the orbitals of each species, 

■ ' ' V'?.^ =0, k',q' ^1,...,Me 

= i{Xk"\Xq") =0, k",q" = l,...,Mc, (47) 

k"q" 

where the differential conditions with respect to the A-spices orbitals have been introduced in 
(25). This leads to the final result for the equations-of-motion for the expansion coefficients: 

"""'(t) = Cf\t) + Cf\t) + Cf\t) + 

+cf''''' (t) + cf^'"'' (t) + cf""^'' it) + cf-^""^ it) + 

+Cf it) + Cf"' it) + Cf"'' it) . (48) 

We remark that other forms of the differential conditions (25,47) can be used, in particular, 
each species can have a different form depending on the physical problem at hand and on 
numerical needs. 

Let us move to the equations-of-motion for the orbitals {0fe(x, i)}, 
{■0fc/(y, t)} and {xfc//(z,i)}. For this, we express the expectation value 



^(ABC)^ in a form which explicitly depends on the various 
integrals with respect to the orbitals. The result is lengthly and posted in Appendix C. 
In particular, the expectation values of the various density operators in H^^^^^ [Eq. (40)] 
emerge as matrix elements of the different intra-species and inter-species reduced density 
matrices. For ease of reading and for completeness, we collect in Appendix B all reduced 
density matrices and their respective matrix elements needed in the theory and its numerical 
implementation. 

We can now proceed and perform the variation of the action functional (45) with re- 
spect to the orbitals. Performing the variation with respect to {(/)^(x, t)}, {V'fc'(y)^)} 
{x^„(z, i)}, making use of the orthonormality relations of each species' orbitals, we solve for 
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the Lagrange multipliers, k,j — l,..., Ma, k' ,j' — 1, . . . , Mb and k" , j" — 1, . . . , Mc'. 

QiA) 



Ma 

E ( 

Ml 

El"!? 

g'=l 



dt 



(^) 



4>q I > 



(49) 



q' / ) 



g"=i 



''dt 



q" 



Xq" 



The terms appearing in the Lagrange multipliers arc all defined in Appendix C. We discuss 
them below, after we arrive at the final form of the equations-of-motion for the orbitals. 
To proceed we introduce the projection operators for the mixture: 

Mb Mc 
u'=l u"=l 

where the projection operator for the 74-species orbitals P*^"^^ was defined in (30). Now, 
eliminating the Lagrange multipliers (49) and making use of the differential conditions for 
each species (25,47), we obtain the final form of the equations-of-motion of the orbitals of 
the mixture, j = 1, . . . , M4, / = !,..., Mb and /' = 1, . . . , Mc: 

Ma , N 



0i 

« \X3") 



p(^) 



(51) 



p(C) 



fc,g=l 

Mb 

fe',9'=l 

Mc 

h^^^ \Xf') + E {P^''\t)}f'k" ( + {P3t^}£J. ) \Xq") 

k",q"=l 



We see the appeahng structure of the equations-of-motion for the orbitals. The various one- 
body operators which assemble the contributions from different orders of the interactions, 
corresponding to the one-, two- and three-body parts of the many-particle Hamiltonian 
^{ABC) (35)^ are separated. Moreover, it is seen that each one-body operator is comprised of 
products of reduced density matrices of increasing order times one-body potentials resulting 
from interactions of the same order (see Appendix C for the explicit terms). This separation, 
originally put forward for the first time in this context for the single-species static theory 
for bosons MCHB [56], is not only theoretically appealing, but is expected to make the 
implementation of the theory in case of higher-body forces further efficient. 
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Equations-of-motion (51) for the orbitals together with (48) for the expansion coefficients 
constitute the propagation theory for mixtures of three kinds of identical particles, interact- 
ing with all possible interactions up to three-body forces. All four possible mixtures (Fermi- 
Fermi-Fermi, Bose-Fermi-Fermi, Bose-Bose-Fermi and Bose-Bose-Bosc) are presented in 
a unified manner, the respective acronyms are denoted as MCTDH-FFF, MCTDH-BFF, 
MCTDH-BBF and MCTDH-BBB. 

To conclude our work, we note that one can compute with imaginary time propagation 
for time-independent Hamiltonians self-consistent ground and excited states for 3-species 
mixtures. Substituting t — > —it into the equations-of-motion for the coefficients and orbitals, 
Eqs. (48,51), the final time-independent (static) theory reads, k = 1, . . . , Ma, k' = 1, . . . , Mb 
and k" = l,...,Mc: 

Ma Ma 

q=l j = l 

Mb Mb 

g'=l j'=l 

Mc Mc 

E Wy^^^'^ + {p^^}i''l" + {pM''1"] m = E /^S" 1^/') ' 

g"=l j"=l 

= £(^^^)Cj-, VJ, (52) 

where, making use of the normalization of the static many-particle wave- function |^(^^^)^, 
^(ABC) _ Y^jCyCj'^^'^^^ is the eigen-energy of the system. Finally, utihzing the fact that the 
matrices of Lagrange multipliers {/i^^^}, {/^iy} and are Hermitian (for stationary 

states) and of the invariance property of the multiconfigurational wave-function (to unitary 
transformations of each species' orbitals compensated by the 'reverse' transformations of 
the coefficients), one can transform Eq. (52) to a representation where {Atfey^}, {/^i^^ 
{a*!"!"} diagonal matrices. This concludes our derivations. 



IV. BRIEF SUMMARY AND OUTLOOK 



In the present work we have specified the MCTDH method for a new complicated system 
of relevance. We have considered mixtures of three kinds of identical particles interacting via 
all combinations of two- and three-body forces. We have derived the equations-of-motion for 
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the expansion coefficients, {Cf{t)}, and the orbitals, {^^(x, t)}, {il^k'{y,t)} and {xfe"(z,i)}, 
see Eqs. (48,51). The self-consistent static theory has been derived as well, see Eq. (52). 

All quantities needed for the implementation of the theory have been prescribed in details. 
On the methodological level, we have represented the coefficients' part of the equations-of- 
motion in a compact recursive form in terms of one-body density operators only, Ipfc^"*!, 
I^iv} I'^i"?"}- "^^^ recursion utilizes the recently proposed Combinadic-based map- 
ping for fermionic and bosonic operators in Fock space [60] that has been successfully applied 
and implemented within the MCTDHB package [62]. Our derivation sheds new fight on the 
representation of the coefficients' part in MCTDHF and MCTDHB without resorting to the 
matrix elements of the many-body Hamiltonian with respect to the time-dependent configu- 
rations, and suggests a recipe for efficient implementation of MCTDH-FFF, MCTDH-BFF, 
MCTDH-BBF and MCTDH-BBB which is well-suitable for parallel implementation. 

As an outlook of the present theory, let us imagine the possibility of conversion between 
the distinct particles, say the conversion of the A and B species to the C species, which can 
be written symbolically as the following "reaction" : 

A + B±^C. 

Such a process would be a model, e.g., for the resonant association of hetero- nuclear ultra- 
cold molecules. The derivation of an efficient MCTDR- conversion theory in this case would 
require the extension of the Combinadic-based mapping [60] to systems with particle con- 
version, and the assembly of more building bricks than just the one-body density operators 
used in the present theory, jp^^^j, jp^yj and 
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Appendix A: Calculating expectation values of operators in mixtures of three kinds 
of identical particles 



Following [60], we write the general expectation value of an operator 0(^"**=^) in a 3-species 
mixture as follows: 



(Al) 



where 

QiSrm.) |^(ABC)(^)^ ^ QiSrm.) ^ | /. = Cf"'^^' {t) \j; t) . (A2) 

J J 

can be a one-, two- or three-body operator or any combination thereof. 
The operation of single-species operators, whether 0^^\ 0*^^^ or 0^*^^ can be read of 
directly from Eqs. (11-19) and we will not repeat them here (one needs just to replace 
therein J a by J in the overall notation, and Ma by Mb or Mq, when appropriate; also see 



For the inter-species two-body operators we prescribe the compact result for 
completeness. For the two-body operators 0(^^) = Sa; fc' « «' ^ifew/^ig^piv' 



Q{AC)_^ q{au) ^[A)au) ^^aq(bc)_y' r,[i^^-) we find- 

^ ~ Z^k,k",q,q" ^kk"qq"Pkq Pk"q" ^ — Z^k',k",q',q" ^k'k"q'q" Pk'q' Pk"q" 

Ma,Mb 



)(AS) ^Pkq Pk'q' 



(i), 



k,k' ,q,q'=l 
Ma,Mc -(a) -(C) 

k,k",q,q"=l 



(A3) 



k',k",q',q"=l 

Note the factorization of the one-body (basic) density operators for the inter-species opera- 
tors, which simplify the way how the coefficients' vector is evaluated. 

For the inter-species three-body operators resulting from the force between two identical 
particles and a third distinct one we list the final result for completeness. For the three-body 
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operators 

A(AAB) _ 1 QiAAB) JA) JB) 

^ ~ 2 ^ ^kk'sqq'lPkslqPk'q'^ 

k,k' ,s,q,q'l 

QiABB) _ 1 V- ri^ABB) 4A)JB) 

^ — 2 2^ ^kk's'qq'l'Pkq Pk's'l'q'^ 

k,k',s',q,q',l' 

QiAAC) _ 1 V O^^^^^ O^^^ cf^^ 
^ ~ 2 2-^ ^kk>'sqq>'WkslqPk"q"-' 

k,k",s,q,q"l 

q(ACC) _ 1 V- rt^ACC) ^{A)JC) 

^ — 2 2-^ ^kk"s"qq"l"Pkq Pk"s"l"q"i 

k,k".s",q,q",l" 

A{BBC) _ I ^(SSC) JB) 4C) 

^ — 2 2-^ ^k'k"s'q'q"l'Pk's'l'q'Pk"q"i 

k'.k",s',q',q"l' 

Q(BCC) _ i q{BCC) ^im^iC) 

^ ~ 2 2-^ ^k'k"s"q'q"l"Pk'q'Pk"s"l"q 

k',k",s",q',q",l" 
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we find: 



. Ma,Mb 



E^iAAB) 
'^kk'sqq'l 



k,k' ,s,q,q' ,1=1 



JA)JB) 

±5siCy ''''' 



4A) 
4A)JB) Pkl 



MaMb 

E^(ABB) 
^kk's'qq'V 

k,k' ,s' ,q,q' ,l'=l 



^(A);,(B) 



±6s>i'C' 



-(B) 

4A)JB) Pyi, 



iP'kq' Pk'q' ^^-J _j_ Qpkq P a' q' 



it) 



. MaMc 



/ , '^kk"sqq"l 



k,k" ,s,q,q" ,1=1 



-yPkq Py'q" . 



4A) 
.(A)4C) Pfci 



c 



q(ACC) 



it) 



Ma,Mc 



2 ^kk"s"qq"l" 

k,k",s",q,q",l"=l 

. Mb,Mc 

\^ n^BBC) 
9 ^k'k"s'q'q"l' 



AC) 



-(C) 

.(.A) AC) Pyiin 
Pkq Ps"q" 



it) 



k',k",s',q',q",l'=l 



.(S) AC) 
_l_r ri'^k'q'Py'q" 

±0.0/;/ o -. 



a(b)a(C) Pyv 

it) T cj'^'''^""" {t) 



c 



Q(BCC) 



it)- 



Mb,Mc 



_ 1 ri(^cc) 

~ 9 2-^ ^k'k"s"q'q"l" 



k',k",s",q',q",l"=l 



.(B) AC) o^~'p^-' 

±5,n"Cy'^''''"'^" it) T c'i'^''^""" (t) 



-(C) 



(A5) 



We remind that the appearance of the one-body (basic) density operators on two levels 
means that the lower-level multiplication has to be performed first, and the upper-level 
second. 

Finally, for the inter-species three-body operator we give the 
closed-form result for completeness. For the three-body operator 

Q{ABC)_Y^ q{ABC) JA)JB)JC) n . 

^ — Z^k,k',k",q,q',q" ^kk'W'qq'q" Pkq Pk'q'PW'q" ^e UUU. 

Ma,Mb,Mc ^a) ab) .(c) 

O(ABC) _ ni^BC) ^Phq Pyq'Py'q" . 



l+\ _ \^ n^ABC) ^Pk. 

y'') ~ ^kk'k"qq'q"^ J 

k,k',k",q,q',q"=l 



it), 



(A6) 



which concludes our Combinadic-based [60] representation of the equations-of-motion for 
the coefficients in MCTDH for mixtures of 3 kinds of identical particles interacting with 
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up to 3-body forces, and the calculations of all relevant matrix elements with respect to 



Appendix B: Reduced density matrices for mixtures of three kinds of identical 
particles interacting with up to three-body forces 

Intra-species reduced density matrices 

The reduced one-body density matrix of the single-species multiconfigurational wave- 
function |*(^n^)> 

is given by: 

p(-^)(xi|x;; t)=NA j 0?X20?X3 • • • d^NA (Bl) 

= *(^)*(x;,X2, . . . , x^^; t)*(^)(xi,X2, . . . ,xjv^; = 

M 
k,q=l 

where its matrix elements in the orbital basis p^^^t) are given in Eq. (28) of the main text. 

Then, the reduced two-body density matrix of the single-species multiconfigurational 
wave- function |*(^^(t)> is given by: 

p(^) (xi , X2 |x; , x'2; t) = Na {Na - 1) y rfx3 ■ ■ ■ d^N^ X (B2) 
X (x'l , x'2, X3, . . . , xjv^ ; t)¥^^ (xi , X2, X3, . . . , xjv^ ; = 

= |{*Ii(x;)*i(x^)*^(x2)*^(xo| = 

M 

= E piliW<^:(x'i,t)</.:(x'2,t)0Kx2,i)0,(Xi,t), 

k,s,l,q=l 

where its matrix elements in the orbital basis p\^}q{t) are given in Eq. (28). 

Finally in the single-species case, the reduced three-body density matrix of \^^'^\t)) is 
given by: 

p(^)(xi, X2, X3|x;, x^, x^; t) = Na{Na - 1){Na - 2) J dx^--- d^N^ x 

x*(^)*(x;, x^, x'3, X4, . . . , xjv^; t)^(^)(xi, X2, X3, X4, . . . , xjv^; t) = (B3) 

= (^(-^)(t) |{*i(x;)*i(x'2)*!,(x'3)*^(x3)*A(x2)*^(xO| *(^)(^))} = 

M 

= E tm^'„ t)<p;{^'„ t)M^s, t)M^2, t)M^^,t), 

k,s,p,r,l,q=l 
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where its matrix elements in the orbital basis p^spriqi^) given in Eq. (28). 

The reduced density matrices of the B and C species are defined in an analogous manner, 
where B and C quantities are to replace the A quantities in Eqs. (B1-B3). 



Inter-species reduced two-body density matrices 

For completeness, we give all inter-species reduced density matrices that occur in a mix- 
ture of three kinds of identical particles interacting with up to three-body forces, where each 
species may have a different spin. There are three such reduced density matrices which are 
associated with the two-body interactions of two distinct particles. 



p^^^) (xi , y 1 1 x'l , y[ ■,t)^ NaNb J o?X2 • • • dx^^ dy2--- dy^g dzi 
x^(^^^)*(x;, . . . , xtv^, y'l, . . . , yivg, zi, . . . , zjv^; t) x 
x*(^^^) (xi, . . . , Xiv^, yi, . . . , y^vs, zi, . . . , zjvc; ^) = 



• dzNc x 



Ma,Mb 



k,k',q,q'=l 



where its matrix elements in the orbital basis are give by: 

.(A) .(B) 



(B4) 



(B5) 



p(^'^)(xi, zi|x'i, z'l; t) = NaNc J (ix2 • • • rfx^v^c/yi • • • dyjVsC/z2 ■ 
^^(ABC)*(^^/^^ . . . , Xat^, yi, . . . , y^^, z'l, . . . , z^^; t) X 
x^'(^^'^')(xi, . . . ,XAf^,yi, . . . ,y7v^,zi, . . .,ZNc;t) = 

= |{*t^(x;)*^(xo*LK)*c(zi)| = 

MaMb 
k,k",q,q"=l 

where its matrix elements in the orbital basis are given by: 



•dZN^ X 



(B6) 



(B7) 
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p'"^^'' (yi, zi|y'i, zi; t) = NbNc J dxi--- dxN^dy2 ■ ■ ■ dyNBdz2 ■ ■ ■ dzNc x 
x*^^^^^*(xi, . . . , xjv^, y'l, . . . , Yat^, z'l, ...,ZNc]t)x 

X^(^^'^)(xi, . . . ,XAr^,yi, . . . ,yArB,Zi, . . .,ZNc;t) = 
Ms, Ms 

= E /^iw WV'fe'lyi- ^)V'.'(yi, ^)x^(z'i, i)x,"(zi, i), (B8) 

k',k",q',q"=l 

where its matrix elements in the orbital basis are given by: 

pii%'At) - E C}^cf'''^'''" it)- (B9) 
/ 

Inter-species reduced three-body density matrices 

There are six reduced three-body density matrices which are associated with the three- 
body interactions of two identical particles with a third distinct one. 



p(^^^)(xi,X2,yi|x;,x'2,y;;t)= (BIO) 
= Na{Na - 1)Nb J dx3--- dxN^dy2 ■ ■ ■ dyNgdzi ■ ■ ■ dz^c x 
x^(^^^)*(x;, X2, . . . , xat^, y'l, y2, . . . , YiVg, Zi, Za, . . . , znc, t) x 
x^'(^^'^)(xi,X2, . . . ,xjv4,yi,y2, . . . ,yArs,Zi,Z2, . . .,ZN^;t) = 

= {¥^^^\t) |{*;,(x;)*^(x'2)*^(x2)*^(xi)*kyi)*B(yi)| = 

Ma,Mb 
k,k',s,l,q,q'=l 

where 

.(A) ,(B) 

are its matrix elements in the orbital basis. 
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where 



p(^^^Hxi,yi,y2|x;,yi,y^;i)= (B12) 
= NaNb{Nb - 1) J dx.2--- (ixjv^c/ys • • • dyNgdzi ■ ■ ■ dz^c x 
x*^^^^^*(xi, X2, . . . , xjv^, y'l, y2, . . . , y^Vs, Zi, Z2, . . . , Zn^; t) x 
^^{ABC) (xi, X2, . . . , xjv^, yi, y2, . . . , yjvs, zi, Z2, • • • , zat^; t) = 

= |{*i(x;)*^(xO*t(y'i)*Uy2)*B(y2)*B(yi)|*(^^^Hi))} = 

Ma, Ms 

= E pilvf4'(^)'^fe(^'i' ^)V'fe'(yi> ^)^My2, (y2, t)V'.'(yi, i), 

fe,fe',s',Z',q',q''=l 

-(A) .(B) 

P&'(0 = E^/W<"'""'"W (B13) 
J 

are its matrix elements in the orbital basis. 



p('^^^)(xi,X2,Zi|x;,x^,z;;t)= (B14) 
= Na{Na - l)Nc J rfx3 • • • dxNAdyi ■ ■ ■ dyNBdz2 ■ ■ ■ dz^c x 
x^(^s^)*(x;, x'2, . . . , xjva, yi, y2, . . . , yjvB, z'l, Z2, . . . , zjv^; t) x 
x^(^^C')(xi,X2, . . . ,XArA,yi,y2, . . . , yAr^, Zi, Z2, . . .,ZN^;t) = 

= {¥^^^\t) |{*i(x;)*i(x^)*^(x2)*4xi)*t(z'l)*c(zi)| ^(^^^)(^))} = 

Ma,Mc 

Yl t)(p*{x.'2, t)(j)l{x2, t)(/'g(xi, t)Xfc»(z'i, t)Xg"izi, t), 

k,k" ,s,l,q,q"=l 

where 

J 

are its matrix elements in the orbital basis. 
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p(^^^)(xi,Zi,Z2|x;,z;,z'2;i) = (B16) 
= NaNc{Nc - 1) J dx2 • • • dxNAdyi ■ ■ ■ dy-Ngdzs ■ ■ ■ dzNc x 
X*^^^^^*(xi, X2, . . . , Xat^, yi, y2, . . . , yjvs, z'l, 4, • • • > ^Nc^t) X 

x*(^^^)(xi, X2, . . . , Xiv^, yi, y2, . . . , y^v^, zi, Z2, ■ ■ ■ , Ziv^; ^) = 

= |{*t^(x;)*^(xi)*t(z;)*t,(zg*c(z2)*c(zi)|*('^^^)w)} = 

pltS'g9"(^)</'fe(xi' ^)Xfc"(z'i, t)x^(z2, ^)X/"(Z2, t)x5"(zi, t), 

fc,fe",s",i",9,9"=l 

where 

= E (0 (B17) 
J 

are its matrix elements in the orbital basis. 



p(^^^nyi,y2,zi|y;,y^,z;;t)= (bis) 

= Nb{Nb - l)Nc J dxi • • • (ixjv^ dys • • • dyNBdz2 ■ ■ ■ dzNc x 
x*^^-^^^*(xi, X2, . . . , x;v^, y'l, y^, . . . , yjvs, z'^, Z2, . . . , z^v^; ^) x 
x^(^^C')(xi,X2, . . . ,XAr^,yi,y2, . . . , y^v^ , Zi, Z2, . . .,ZN^;t) = 



Mb,Mc 

E Pm'lq'q"itWk'iy'l, tWs'iy2, 0^i'(y2, ^)</'g'(yi, ^)Xfc"(z'i, t)Xg"{zi, t), 

k',k",s',l',q',q"=l 

where 

pil^'sl^A^) - E ^X^)^^?'"'^""'" (^) (B19) 
J 

are its matrix elements in the orbital basis. 
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p(''^''nyi,Zi,Z2|y'i,z;,z^;t)= (B20) 
= NbNc{Nc - 1) J (ixi • • • dxN^dy2 ■ ■ ■ dyNgdzs ■ ■ ■ dz^c x 
x*^^^^^*(xi, X2, . . . , xjv^, y'l, y2, . . . , yNB,z[, z'^, . . . , zn^; t) x 
x*(^^^)(xi, X2, . . . , xjv^, yi, y2, . . . , y^v^, zi, Z2, • • • , zat^; ^) = 

= |{*Uyi)*B(yi)*k<)*kz2)*c(z2)*c(Zl)|*(^^^)W)} = 

XI Plwr'9'g''(*)V'fe'(y'l, ^)^g'(yi, ^)Xfc"(z'l, ^)X^(Z2, ^)X/"(Z2, ^)Xg"(Zl, t), 

k',k",s",l",q',q"=l 

where 

.(S) -(C) 

piS".'."(t) = (B21) 
/ 

are its matrix elements in the orbital basis. 

Finally, there is a single reduced three-body density matrix which is associated with the 
three-body interaction of three distinct particles. 



p(^^^)(xi,yi,zi|x;,yl,z;;t)= (B22) 
= NaNbNc J dx2 • • • dxNAdy2 ■ ■ ■ dyNBdz2 ■ ■ ■ dzNc x 
x^(^^C')*(x;,x2, . . . ,xjv^,yi,y2, . . . , y^v^, z'^, Z2, . . • ,zjvc;t) x 

X^(^^'^)(Xi,X2, . . . ,XAr^,yi,y2, . . . ,yiv^,Zi,Z2, . . .:ZNc;t) = 

= {¥^^^\t) |{*i(x;)*^(xO*kyl)*B(yi)*t(4)*c(zi)| = 

Ma,Mb,Mc 
k,k',k",q,q',q"=l 

-(A) -(C) 

pii^^'L'it) = E c}{t)c';^ (t) (B23) 

are its matrix elements in the orbital basis. 



where 



33 



Appendix C: Further details of the derivation of the equations-of-motion for mix- 
tures of three kinds of identical particles 

The derivation of the equations-of-motion for the orbitals (51) starts from expressing the 
expectation value of H(^^'^^ with respect to the many-particle wave-function |^(^^'^)^ in a 
form which depends explicitly on the various integrals with respect to the orbitals. Thus we 
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have: 



' k,q=l 



dt 



kq 



+ 



-.Ma ..Ma 
+ 1 V n^^^W^^^+- T 0^^^ U^^^ + 

' ^ / ^ rkslq ksql ' g / j rksprlq kspqlr ' 

fe,s,p,r-,i,q=l 



+ 



k,s,l,q=l 

Mb 

fc',9'=l 



^'^v Ydt 



Q(By 



+ 



-.Mb Ms 

+ 1 V o^^^ W^""^ +- V o^^^ c/f''^ + 

^2 f^k's'l'q'^^ k's'q'V ^ Q /-^ Pk's'p'r'l'q''^ k's'p'q'l'r' ^ 

k',s',l',q'=l k',s',p',r',l',q'=l 

r , 1 



+ E Pi 



k",q"=l 



"'k"q" S ' 



at 



+ 



+ 



1 



Mc 

/ , Pk"s"l"q"*^ k"s"q"l' 
:",l",q"=l 



k"q"_ 
+ 



2 ^ 

k",s",l",q"=l 

-t- Pk"s"p"r"l"q"'^k"s"p"q"l"r" ^ 

k",s",p",r",l",q"=l 

Ma, Mb MaMc 
V- (AB) ^(AB) (AC) yr^(AC) 

^ Pkk'qq'^^ kk'qq' ^ Pkk"qq"^^ kk>'qq" ^ 

k,k',q,q'=l k,k",q,q"=l 

Mb,Mc 

^ Pk'k"q'q"^^ k'k"q'q" ^ 

k',k",q',q"=l 

Ma,Mb . Ma, Ms 

-ui ^(^^5) 7-r(^^B) , i ^(^SS) ..(ABB) 

"•"2 Pkk'slqq''^kk'sqq'l^ 2 Z^ P kk' s' V qq' kk' s' qq' V ^ 

k,k',s,q,q',l=l k,k' ,s' ,q,q' ,l'=l 

Pkk"s"l"qq 



1 



Ma.Mc Ma,Mc 

/ V Pkk"slqq"'-'kk"sqq"l 2 Z^ 
i:",s,g,g",J=l k,k",s",q,q",l"=\ 



Ma,Mc 

/ , Pkk" a"!." nn"^ kk" s"qq"l" 



MbMc 

/ V Pk'k"s"l"q'q"'^k'k"s"q'q"l" 

■ c/t „l f^ll 111 — 1 



"2 „i:„ 

A;,/c ,s,q,q ,1 — 1 5^ ,v 

Mb,Mc Mb,Mc 

^2 Z^ Pk'k"s'l'q'q"'^k'k"s'q'q"l' ^ ^ Z^ 

k',k",s',q',q",l'=l k' ,k" ,s" ,q' ,q" ,1" = 

Ma,Mb,Mc 

MBC) ..(ABC) _ 
^ Z^ Pkk'k"qq'q"'^kk'k"qq'q" / , ''^ jV' )^ JV' ) ■ 

k,k' ,k" ,q,q' ,q"=l |j| 

density operators appearing 



(CI) 



{J} 

in (CI) have been pre- 



The expectation values of the various 
scribed in Appendix B. 

The matrix elements in (CI) of the A and correspondingly of the B and C single-species 
terms with respect to the orbitals have been discussed in section II A, see Eq. (4). The 
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matrix elements arising from two-body inter-species interactions are listed for completeness 
below: 

W^i'W = // rAy,t)xU^,t)W^''^\y,z)^p,,{y,t)x,''{z,t)dydz, (C2) 
and the matrix elements arising from three-body inter-species interactions read as follows: 



rj{AAB) ^ 

kk'sqq'l 



(/**(x, t)0:(x', t)rAy, t)m^^^\^, x', y)0,(x, t)M^', t)^,Ky, t)d^d^'dy, 



jAABB) 
'^kk's'qq'l' 




jAAAC) 
^kk»sqq"l 




(x, t)rAy, t)rAy', t)^(^^^)(x, y, y')0,(x, t)^,,(y, t)^,(y', t)d^dydy', 



:(x, t)0:(x', t)x^(z, i);7(^^^)(x, x', z),^,(x, t)0,(x', i)x,"(z, t)dxdx'(iz, 



tAACC) 

'-^kk"s"qq"l" 



0^(x, f)x^(z, i)x:"(z', t)^(^^^)(x, z, z')<^,(x, i)x5"(z, i)x/"(z', t)d^dzdz', 



jABBC) 
'^k'k"s'q'q"l' 



Ak'iy, t)rAy', t)xl"{^, t)u^''''''\y, y', z)My, t)My', t)xq"{^, t)dydy'dz, 



jABCC) 
'^k'k"s"q'q"l" 



^ III ^fe'(y'^)^^(^'^)^>(^^^)^^''''''^(y'^'^0V'.Ky,^)x5'Kz,^)xp'(z',t)dycizdz', 



rAABC) _ 

^kk'k"qq'q" ~~ 



^ III '^fe(^'^)^fc'(y'^)^^(^'^)f^^^''''Hx,y,z)(/',(x,t)^,,(y,t)x,"^,^)rfxdydz. 

(C3) 

Performing the variation of the integrals (C2) with respect to the orbitals {0jk(x, t)}, 
{ipk'{y,t)} and {xa;"(z,^)}, we find six types of inter-species one-body potentials emerging 
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from two-body interactions: 

wS^\y,t) = 1 0Ux,Oiy(-^^)(x,y)0,(x,t)(ix, 

^t'f\^,t) = J rAy,t)W^''''\y,z)ij,,{y,t)dy. (C4) 

Making the variation of the integrals (C3) with respect to the orbitals, we arrive at fifteen 
types of inter-species one-body potentials resulting from three-body interactions: 
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uit'tfM- J I 0:(x',t)#(y,t)f/(^^^)(x,x',y)0,(x',i)V;,'(y,^)^x'(iy, 
Uii^^\y,t) = 1 1 0^(x,t)</.:(x',i);7(^^^)(x,x',y)0,(x,t)(/,Kx',t)rfxdx', 

ujff^jM^ J I rAyM'{y\t)u^''^^\^.y,y')MyMi'{y'.t)dydy\ 

uiXP{y,t)^ J 1 0:(x,i)^:,(y',i)t/(^^^)(x,y,y')0,(x,i)V;ay',t)cixdy', 
t>iw(x,t) = 1 1 0:(x',t)x^(z,t)t/(^^^)(x,x',z)(/),(x',Ox,"(z,t)dx'(iz, 
= II </>*(x,t)C(x',t)[/(^^^'(x,x',z)(/),(x,t)(/)Kx',t)c^xdx', 

Ul%fj{z,t) = 1 1 (/»^(x,t)x:.(z',t)t/(-"^'^)(x,z,z')0,(x,t)xr(z',t)dxdz', 

^iSv'(y,^) = / / ^My',Ox^(z,t)f/^^^^Hy,y',z)^,(y',^)x."(z,Ociy'ciz, 
= / 1 #(y,t)^My',t)c/(^^^ny,y',z)V'.'(y,^)V'^'(y',i)c^yc^y', 

t>wV'(y'^) = / / X^(z,^)x:"(z',i)^^^^^ny,z,z')x,"(z,^)XF'(z',^)c?zdz', 

C"5''(z,0 = / 1 ^:^(y,t)x:"(z',t)^7(^^^ny,z,z')^,'(y,^)x/"(z,0c?yc^z', 

= // V'fc'(y,t)x^(z,^)C>^^^''nx,y,z)z^,Ky,i)x,"(z,t)c/yc/z, 
t^w(y,^) = // 0:(x,t)x^(z,i)f7(^^^nx,y,z)<^,(x,t)x,"(z,t)cixciz, 

= // </':(x,t)#(y,t)f7(^^^)(x,y,z)(^,(x,0V','(y,0^^xdy. (C5) 

All one-body potentials in (C4) and (C5) are local (for spin- independent interactions), time- 
dependent potentials. 

To arrive at the final form of the equations-of-motion (51), we define the auxiliary one- 
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body operators for the A-species' particles: 

Ma Mb Mq 

iP-i^ hq = Z^Pkslq^sl + 2^ Pkk'qq'^k'q' + 2^ P kk" qq" ^ k" q" ^ 

s,l=l k',q'=l k",q"=l 

Ma Ma,Mb 

in /7\(^) = i n(^) /7(^) ^ rl-^^^^ jMAB) 

\P3'-^ fkq "2 Pksprlq^ splr ' 2-^ Pkk' slqq'^ sk'lq' ' 

s,p,l,r=l k',s,q',l=l 

Mb MaMc 

Sr^ (ABB) jAABB) sr^ (AAC) fAAAC) 

2-^ Pkk's'l'qq''^k's'q'l' ^ Pkk"slqq"'^ sk"lq" 

k',s',q',l'=l k",s,q",l=l 

Mc MbMc 

, MCC) jAACC) (ABC) jAABC) /p^x 

"I" Pkk"s"l"qq"^k"s"q"l" Pkk'k"qq'q"'^ k'k"q'q" ' \^^/ 

k",s",q",l"=l k',k",q',q"=l 

for the i?-species' particles: 

Mb Ma Mc 

TT,T.(B) _ (B) yiAB) , \^ (AB) yrABA) (SC) t9.(BC) 

s',l'=l k,q=l k",q"=l 

Mb Ma 

XP^'^ ik'q' — 2 Pk's'p'r'l'q''-' s'p'l'r' ^ Pkk'slqq'^ ksql 

s' ,p' ,r' ,1' =1 k,s,q,l=l 

Ma, Mb MbMc 

Pkk's'l'qqi'^ks'ql' Z^ Pk'k"s'l'q'q"'^ s'k"l'q" 

fe,s',g,i'=l fc",s',9",;'=l 

, IBCC) jABCC) (ABC) .7(BAC) 

Z^ Pk'k"s"l"q'q"'^k"s"q"l" ^ Z^ Pkk'k"qq'q"'^ kk"qq" ^ \^ ' ) 

k",s",q",l"=l k,k",q,q"=l 

and for the C-species' particles: 

Mc Ma Mb 

\P'2^^ Jk"q" — Z^ Pk"s"l"q"^^ s"l" ^ Pkk"qq"^^ kq Z^ Pk'k"q'q" k'q' ' 

s",i"=l fc,q=l jfc',q'=l 

Mc Ma 
lA'SL' — 2 ^ Pf.„g„p„j.,q„g„Ug„p„i„^„ -\- Pkk" slqq"^ ksql 

s",p",r",i"=l fe,s,g,i=l 

MaMc Mb 

E iACC) jACAC) ^(BBC) frCCBB) 

Pkk"s"l"qq"'^ks"ql" ^ Z^ Pk'k" s'l'q'q"'^ k' s'q'l' 
k,s",q,l"=l k',s',q',l'=l 

Mb,Mc Ma,Mb 

E iBCC) jACBC) Sr^ (ABC) fACAB) 

Pk'k"s"l"q'q"'^k's"q'l" Z^ P kk' k" qq' q" kk' qq' ■ \^°) 

k',s",q',l"=l k,k',q,q'=l 

These auxiliary one-body operators are constructed from products of matrix elements of 
reduced density matrices of increasing order (see Appendix B) times the one-body potentials 
resulting from interactions of the same order, see Eqs. (C4) and (C5). The derivation of the 
equations-of-motion (51) is now fully completed. 
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